Temperature dependent Bogoliubov approximation in the classical fields approach to 

weakly interacting Bose gas 



o 
o 

(N 
> 

o 

(N 



> 

(N 
(N 



o 

c3 



i 

c 

o 
o 



Miroslaw Brewczyk, 1 Peter Borowski, 2 Mariusz Gajda, 3 and Kazimierz Rzazewski 2 
1 Uniwersytet w Bialymstoku, ulica Lipowa 1^1, 15-424 Bialystok, Poland 
2 Centrum Fizyki Teoretycznej PAN, Aleja Lotnikow 32/46, 02-668 Warsaw, Poland 
3 Instytut Fizyki PAN, Aleja Lotnikow 32/46, 02-668 Warsaw, Poland 
(Dated: February 2, 2008) 

A classical fields approximation to the finite temperature microcanonical thermodynamics of 
weakly interacting Bose gas is applied to the idealized case of atoms confined in a box with periodic 
boundary conditions. We analyze in some detail the microcanonical temperature in the model. We 
also analyze the spectral properties of classical amplitudes of the plane waves - the eigenmodes of 
the time averaged one-particle density matrix. Looking at the zero momentum component - the 
order parameter of the condensate, we obtain the nonperturbative results for the chemical potential. 
Analogous analysis of the other modes yields nonperturbative temperature dependent Bogoliubov 
frequencies and their damping rates. Damping rates are linear functions of momenta in the phonon 
range and show more complex behavior for the particle sector. Where available, we make comparison 
with the analytic estimates of these quantities. 



I. INTRODUCTION 

Since the achievement of Bosc-Einstcin condensation 
in dilute atomic gases there is a remarkable experi- 
mental activity in this area of physics. Quantum proper- 
ties of the weakly interacting Bose gas are studied both 
close and away from thermal equilibrium. While in most 
cases the experiments are performed with as cold gas 
sample as possible, some experiments test the depen- 
dence of these properties on temperature. The central 
role in the theoretical studies of this system is played by 
the collective excitations. At and near zero temperature 
these collective excitations are described very well by the 
Bogoliubov approximation. Much less is known about 
the collective excitations at higher temperatures. 

More generally, quantum theory of weakly interacting 
Bose gas at finite temperatures is still a challenge. To 
answer this challenge most researchers describe the sys- 
tem as consisting of two different, mutually interacting 
components: the condensate and the thermal cloud |2|. 
This way a significant progress has been achieved in ex- 
plaining such experimental observations as, for instance, 
temperature shifts and damping rates of various oscilla- 
tion modes of the two component system The two 
component approaches, although successful, are unsatis- 
fying in their underlying assumptions. The splitting of 
the gas at nonzero temperature into the condensate and 
the thermal cloud should be the result of the Bose statis- 
tics and of interactions. 

In a series of papers several groups formulated the, so 
called, classical fields approximation jj, |jj |(J 0, H, H, [H| • 
Within this formulation there are, however, two slightly 
different approaches. The group of ENS developed 
the so called truncated Wigner method for Bose con- 
densed gases. The main idea is to describe the total 
system by a classical field obeying the Gross-Pitaevskii 
equation. The classical field, except of a condensate com- 
ponent, contains also a stochastic part representing a 
thermal cloud. Mean values of observables are calculated 



as averages over an ensemble of classical fields which are 
chosen according to the Wigner quasi-distribution func- 
tion of the initial thermal equilibrium density operator of 
the gas. This approach corresponds to the canonical de- 
scription of the system where a temperature not energy 
is a control parameter. 

The approach that we use has been formulated in 0, 
and corresponds to the microcanonical description. This 
approximation, drawing from quantum optics, consists 
of replacing the quantum mechanical operators of highly 
occupied modes of the Bose field with complex c-number 
amplitudes. The idea of replacing the annihilation (cre- 
ation) operators by c-number amplitudes is a straight- 
forward generalization of Bogoliubov hypothesis to a case 
when large number of modes is macroscopically occupied. 
It is the aim of this paper to analyze the classical fields 
approximation in terms of such fundamental notions as 
the microcanonical temperature, the chemical potential 
and the quasiparticle excitations. 

The paper is organized as follows: In Section [HJ for 
the sake of completeness, we briefly describe the classical 
fields approximation. In Section 1 1 1 1 1 we analyze the dy- 
namical equations of the approximation from the point of 
view of underlying spectra. In Section llVI the numerical 
results are discussed to support the analysis of Section 
IIIII In particular we show that the Bogoliubov energies 
at high temperature follow the gapless formula derived in 
[TT| . We also analyze the thermal damping rates of the 
modes. We find that they are proportional to momenta 
for the phonon excitations. The dependence of damping 
rates of particle-like modes on momentum is more com- 
plex showing the existence of bending points. We present 
some concluding remarks in Section Ivl 



II. CLASSICAL FIELDS APPROXIMATION 

We consider N identical atoms subject to Bose- 
Einstein statistics, confined to a cubic box of volume 
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V , with their atomic wave functions satisfying periodic 
boundary conditions. The atoms interact via contact po- 
tential characterized by the scattering length a s . We use 
the method of second quantization. Our dynamical vari- 
able is a bosonic field operator ^(r,t) that destroys a 
particle at position r and obeys standard bosonic com- 
mutation [*(r,t), ^(r'jt)] = 8(r - r') relations. The 
Hamiltonian of the system reads: 



H = f d 3 r ¥(r,t)—V(r,t) 
J 2m 

■ J d 3 r * f (r, i)* f (r, t)#(r, t)V(r, t) . 



(1) 



The Heisenberg equation of motion for the field oper- 
ator follows: 



&(r,t)${r,i)${r,t) . (2) 



We do not know how to solve this nonlinear operator 
equation. However, a natural simplification is possible 
for the degrees of freedom (modes) of the field, which 
are occupied by large number of atoms. Just like corre- 
sponding highly occupied modes of electromagnetic field, 
they can be described by c-number complex amplitudes 
rather than by their creation and annihilation operators. 
On the other hand the modes that are sparsely popu- 
lated and require full quantum treatment, may be in the 
crudest approximation neglected. 

The symmetry of the box with periodic boundary con- 
dition determines a natural set of modes for the ther- 
mal equilibrium states of atoms. The natural modes 
are simply the plane waves, with quantized momentum 
p = 2irh(ni,Ti2,n3)/L (n, being integer). Therefore, the 
field operator can be expanded in these modes: 



*(r,<) = -^=^exp(-ipr/ft)a p (t) 



(3) 



Substituting the expansion J3J into (J2J we get a set of 
nonlinear equations for the creation and annihilation op- 
erators of the plane wave modes: 



da p (t) 
dt 



2mh 



^ Yl «qi a q2«P+qi-q2 , ( 4 ) 

qi,q2 



where the coupling g is equal to: 

4irh 2 a s 
9 = 



(5) 



At this point we identify the set of highly occupied 
modes. By definition these are all long wave length de- 
grees of freedom up to a maximal cut-off momentum 



Pmax- For these modes we replace the operators with 
complex amplitudes: 



(6) 



in such a way that the square of the modulus of the am- 
plitude gives the fraction of atoms in a given mode. The 
above replacement is a direct generalization of the Bo- 
goliubov hypothesis. 

Confining our attention to the highly occupied modes 
we are reducing the dynamical equations to the set of 
nonlinear differential equations for the amplitudes: 



da p (t) 
dt 



= —i 



2mh ap 



■gN \- * 

^ 2s a qi a q2« P - 

qi,q2 



qi-q2 



(7) 



which may be solved numerically. Note that summa- 
tion over momenta is truncated at |qi| = |qa | = p m ax, 
where p max is a value of momentum of the highest mode 
which occupation is still macroscopic. In numerical im- 
plementation of the method the value of p ma x has to be 
determined self-consistently. 

Eq. is a momentum representation of the stan- 
dard Gross-Pitaevskii equation. For numerical purposes 
it is convenient to switch to position representation which 
corresponds to the c-number version of the operator Eq. 
|j2J on rectangular grid with the grid spacing defined by 
a cut-off momentum Ax = 7r/p max . In the classical 
fields method, a finite grid version of the Gross-Pitaevskii 
equation describes a total system - Bose-Einstein con- 
densate and a thermal cloud. 

This interpretation is shared by other groups which 
apply the classical fields method for microcanonical de- 
scription of a cold Bose system 0, 0, H ■ In the recent 
paper the value of the cut-off momentum has been 
chosen arbitrary. It means that only a fraction of all 
macroscopically occupied modes is taken into account, 
therefore in such a treatment there is no simple method 
of determination of a number of particles in the whole 
system. This leads to difficulties in assigning an abso- 
lute value of temperature and authors of [3] calculate 
some scaled temperature T/N rather than T. We treat 
the problem of a cut-off momentum differently. In the 
present approach a number of classical modes, or equiv- 
alently a cut-off momentum, is an important physical 
parameter. We determine its value by requiring that an 
occupation of the cut-off momentum mode is equal to 
one. This is a kind of compromise as the validity of ap- 
proximation is limited to modes with occupation large 
compared to one, but on the other hand we want to de- 
scribe the whole system. Such a choice of the cut-off 
momentum ensures that the total occupation of modes 
with momenta not included in our computations is only 
a few percent even close to critical temperature. More- 
over, we are able to determine the absolute (not scaled) 
value of temperature of the system. 

The most important finding is 0, Q that almost any 
initial condition of given number of atoms and given en- 
ergy, after the transient thermalization period leads to 
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a steady state that does not depend on the particular 
choice of initial conditions. This way we are able to 
produce a numerical version of fluctuating, weakly in- 
teracting Bose gas in the quantum degenerate region. Of 
course, at temperatures close to the critical one there 
are many highly occupied modes of the Bose field. All 
authors using the classical fields approximation identify 
the zero momentum component of the wave function with 
the condensate and all other components as a part of the 
thermal cloud. In a recent paper [9j we have investigated 
this problem and came to the conclusion that the coarse 
graining due to the finite resolution of the realistic mea- 
surement, both in time and in space, reduces the pure 
state described by a single wave-function of the method 
to a mixed state in which there is no coherence between 
various momentum components of the field. Only with 
this interpretation the perfectly isolated system of a mi- 
crocanonical approach may be viewed in quantum me- 
chanics as a mixed state. In the next Sections we are go- 
ing to spectrally analyze the long time solutions of Eqs. 
(0) to generalize the celebrated Bogoliubov approxima- 
tion to higher temperatures. 



III. EXCITATION SPECTRUM: 
APPROXIMATE TREATMENT 

In this section we present an approximate analytical 
treatment of dynamical equations of the classical fields 
method. Under assumption of a steady-state evolution 
we discover some constants of motion what allows us to 
linearize the equations and find an excitation spectrum of 
the system. This approach leads to the notion of quasi- 
particles and the famous Bogoliubov-like spectrum. 

The equations for complex amplitudes originating from 
the many body Hamiltonian are given by Q) and in a 
particular units of length (the size of the box L) and 
time (mL 2 /h) are: 



l,m 



(8) 



The set of nonlinear Eqs. © has two constants of 
motion. The first is the total energy of the system, 
E = iVE k (k 2 /2)n k + (ffJV/2)Ek,i, m tti«f«m« k+ i-m 
and the second is the number of particles N ^ k rik, where 
n k = a k a k- A s we have already mentioned, the system 
reaches a thermal equilibrium [f| after some transient 
time, which depends on constants of motion only. In this 
state, a condensate occupation undergoes small fluctua- 
tions around a well defined mean value. Neglecting these 
small fluctuations we assume that uq = ctQa is constant. 
A value of rig depends on the energy of the system and 
is close to one for small energies and decreases with in- 
creasing energy. 

By rearranging terms in Eq. JSJ we obtain the follow- 
ing equation for the amplitude of the condensate (p = 0) 



mode: 

«o = ~igN(2 - n )a - igN ( ^ "kO-k 

.k/o 



where 



fo(t) = -igN 2J a k aq"k- 

k^l;k,l#0 



fo(t) , 
(9) 

(10) 



plays a role of an external force. Eqs. 10 and its com- 
plex conjugate describe a pair of coupled harmonic os- 
cillators driven by external forces. The coupling term, 
Sk^o a kC*-k, is an anomalous density. 
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FIG. 1: Fourier transform (in arbitrary units) of the anoma- 
lous density (peaks in the right part of the figure) and driving 
force (peaks in the left part of the figure) for various conden- 
sate occupations as indicated in the legend. 

In general time dependence of both the anomalous den- 
sity and driving force can be very complicated. However, 
in a steady state both these quantities simply oscillate in 
time: 



k^o 

fa = -K^e-^+^ 



(11) 
(12) 



where fi plays a role of a chemical potential (as it can 
be seen later), <f) is some constant phase while 6 and n 
are (real) constants of motion depending on the energy 
and number of particles only. To support the assumed 
ansatz in Fig. ^ we show the Fourier transforms of the 
anomalous density and the driving force for three differ- 
ent values of the total energy. Indeed, independently of 
the energy, the Fourier transform of the anomalous den- 
sity is sharply peaked at frequency which is larger by the 
factor of two than frequency corresponding to a peak in 
the spectrum of the driving force. Because of Eqs. 
and (|12l) the equation for condensate amplitude has a 
simple solution: 



n e 



(13) 
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what is consistent with our assumption that condensate 
population n Q does not depend on time. Evidently, fj, 
plays a role of a single particle energy in a condensate 
phase. The Eq. gives also relation between chemical 
potential [i and other, introduced above parameters: 



(j, = gN(2 - n ) - gN(S + k) 



(14) 



The anomalous density <5 is small at low energies or close 
to a critical energy, having maximum in between. The 
parameter k is very small at low energies but it grows 
continuously and gives a significant correction to the 
chemical potential at large energies. Our result can be 
compared to the one obtained from the two gas model 
fi = gN(2 — no) (see Fig. |2J). Only at very small en- 
ergies this model predicts a value of chemical potential 
correctly. 

In Fig. [21 we show the chemical potential as a func- 
tion of a condensate occupation (which is a monotonic 
function of temperature). We present two curves, both 
corresponding to the same value of the effective coupling 
gN = 731.5, but for two different grid sizes. The number 
of grid points is equal to a number of macroscopically 
occupied modes. If we increase the number of macro- 
scopically occupied modes while keeping population of 
the highest mode constant, we inevitably increase the 
total number of particles in the system, and simultane- 
ously decrease the value of g (gN is fixed) There- 
fore, the curve corresponding to the larger grid (circles) 
describes larger system with weaker interactions as com- 
pared to the curve marked by boxes. In Fig. [3]we present 
the chemical potential as a function of the interaction 
strength g at fixed condensate fraction no = 0.55 and 
fixed effective coupling gN — 731.5. According to the 
two gas model the chemical potential should be constant 
in such a case and its value, for chosen parameters, equals 
/i = 1060. This value is indicated in Fig. Elby horizon- 
tal dashed line. Our results show that chemical potential 
depends on the interaction strength not only through the 
product gN. We expect that two gas model result can be 
reached in the limit of g —> 0, however calculation in this 
regime is a demanding task. In the inset we show the 
chemical potential versus condensate population (tem- 
perature) for fixed value of both: particle number N and 
interaction strength g. The difference between the two 
gas model prediction and our result grows with temper- 
ature. 

Now, we rearrange terms in the dynamical equation 
© for amplitudes of excited modes to get the following 
expression: 



% + 2gN j a p 

igN \al + ^2 a k«-k ] a*_ p + f p , (15) 
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FIG. 2: Chemical potential as a function of the condensate 
population for the effective coupling strength gN = 731.5. 
Circles and boxes show numerical results obtained on the grids 
36x36x36 and 16x16x16, respectively. The solid line comes 
from the formula derived within the two gas model. 
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FIG. 3: Chemical potential as a function of the interaction 
strength g alone for gN — 731.5 and the condensate oc- 
cupation no = 0.55 (solid line is added to guide an eye). 
The horizontal dashed line indicates the value of the chem- 
ical potential obtained from the two gas model formula: 
fj, = gN(2 - n ) ~ 1060. The inset shows, just like Fig. H the 
chemical potential as a function of the condensate population 
but for a constant number of atoms N = 235000. 



where 

f P = -ig N ^2 ak a i"p+k-i (16) 
M-p1;¥p 

is a driving force of the a p mode. Guided by the previ- 
ous experience, we expect that the driving force / p be- 
comes important at larger energies. However, its time 
dependence is not as simple as that of /o. In the first 
approach we omit this force in our analytical treatment. 
We are well aware that doing this we loose very impor- 
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tant contribution to the dynamics, nevertheless even this 
oversimplified analysis gives interesting physical insight. 

After all the above simplification we get the set of 
coupled equations for amplitudes of a p and al p modes. 
They describe excitation of the system by creation of cor- 
related pair of particle with momentum p and a hole in 
the condensate with opposite momentum. The solution 
of this set reminds famous Bogoliubov transformation: 



„„ = e f e-^Pp + e*"* _^ f _ " P ^ Pt p ) , (17) 



gN(n - 5) 



where amplitudes p p and /?* p are some constants, and 
e p is the temperature dependent Bogoliubov spectrum: 



and u>p is: 



e p = ja;* - [gN(n - 8W 



u) p = — + gN(n + k + 6) . 



(18) 



(19) 



The amplitude of the a p mode is a superposition of two 
components: the one oscillating with frequency /i + e p 
and the second with frequency [i — e p . Excitations are 
Bogoliubov quasiparticles. We see, that this well known 
fact is reproduced by the classical fields method and is 
also valid at higher temperatures and for large interaction 
strength (in fact it is valid up to a critical temperature 
what we will show in the next section) . At low momenta 
both frequency components of the amplitude a p are com- 
parable while at high momenta the amplitude of the neg- 
ative frequency component vanishes. Eq. I|18|l predicts a 
gap, i.e., the energy of excitation does not tend to zero 
as the momentum p tends to zero. This result violates 
the Hugenholtz-Pines ^3] theorem which shows that ex- 
citation spectrum is gapless. The reason is that in our 
analysis we have neglected the term / p . Inclusion of this 
term gives the correct gapless spectrum and allows for 
determination of Landau damping rates of excitations. 



IV. NUMERICAL RESULTS 

In this section we show results obtained from numeri- 
cal solution of the dynamical Eqs. ||SJ and compare them 
with approximate analytical solutions. We focus on the 
thermodynamics of a uniform weakly interacting Bose 
gas in the equilibrium. Knowing the spectrum of the ele- 
mentary excitations (quasiparticles) as well as their pop- 
ulations we deliver the scheme which assigns the temper- 
ature to the system. After that most of the thermody- 
namic properties of the system can be derived, for exam- 
ple, the condensate occupation as a function of tempera- 
ture. On the other hand, the width of the quasiparticles 
originating from the finite lifetime of the elementary ex- 
citations allows one to investigate dissipative processes 
in the system. 



As we have already mentioned in Section [D] instead 
of solving Eqs. JSJ it is more convenient to trans- 
form them to the position representation what leads 
to the finite- grid version of the time-dependent Gross- 
Pitaevskii equation for a system of particles in a box with 
periodic boundary conditions. The equation is written 



(20) 



where particular units for the length (which is the box 
size L) and the time (mL 2 /h) are assumed and the cou- 
pling constant g equals Ai:a s /L (a s being the scattering 
length). This equation has to be solved on a rectangu- 
lar grid and the grid spacing Ax determines the value 
of the cut-off momentum p max — it/ Ax. The initial 
wave function is generated from the ground state solu- 
tion by its random disturbance followed by the normal- 
ization. The strength of the disturbance determines the 
total energy per particle. Both the energy and the par- 
ticle number are preserved during the evolution of the 
Gross-Pitaevskii equation and the system reaches the 
same stationary state independently of the choice of the 
initial wave function provided that the energy and the 
particle number are kept constant. 
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FIG. 4: Fourier transform (in arbitrary units) of the (0,0,1) 
and (3,3,3) modes for gN=731.5 and no = 0.79 (low temper- 
ature). Note the existence of two groups of frequencies. 

As it was discussed in the previous sections, the high 
energy solution of the Gross-Pitaevskii equation rep- 
resents a pure state of the system and only time av- 
eraging procedure related to the nature of the obser- 
vation process allows one to split the system to con- 
densed and non-condensed parts. To this end, the time 
averaged (after the system enters the stationary state) 
single-particle density matrix is diagonalized and coher- 
ent modes (among them the condensate mode) and their 
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FIG. 5: Fourier transform (in arbitrary units) of the (0,0,1) 
and (3,3,3) modes for gN=731.5 and no = 0.07 (high temper- 
ature) . 



populations are obtained. For a uniform gas the eigen- 
functions of the single-particle density matrix are just the 
plane waves and the condensate is identified as the mode 
with zero momentum. Therefore, we propagate the high 
energy solution of the Gross-Pitaevskii equation and, at 
each time step, do the decomposition into plane waves by 
using the Fourier transform technique in spatial domain. 
In such a way we gain the time dependence of natural 
modes. Their spectrum (Fourier transform with respect 
to time) is plotted in Figs. 0] and [5] in the case of low and 
high temperatures, respectively, for modes (0,0,1) and 
(3,3,3). Typically, two groups of frequencies are visible 
in the spectrum and the separation between them in- 
creases with the mode energy. Moreover, increasing the 
energy of the mode leads also to the suppression of the 
lower frequencies group. Another important property is 
that the width of each group is getting broader when the 
temperature of the system is increasing. 

The mean frequency of each group of frequencies is well 
described by the gapless Bogoliubov-like formula: 




gNn ) -(gNnoY 



(21) 



where no is the condensate fraction. This simple general- 
ization of the Bogoliubov formula is known as the Popov 
approximation [llj. In Fig. H3 which is the case of (1, 1, 1) 
mode, frequencies [i ± e p are marked by the vertical dot- 
ted lines whereas those given by the simplified expression 
(|18fl (when the driving force term is neglected) derived in 
the previous section by the dashed ones. Obviously, the 
formula i|21[l works better. The excellent agreement be- 
tween the Bogoliubov-like expression l|21|) and numerical 
spectrum manifests itself in Fig. \7\ where we plot the 



quasiparticle energies for various temperatures (conden- 
sate populations) and interaction strengths. Numerical 
points (marked by boxes, circles, and triangles) are cal- 
culated as the mean frequency values (with respect to fi) 
corresponding to higher frequency groups, whereas the 
solid lines come from (|21|) . 
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FIG. 6: Fourier transform (in arbitrary units) of the (1,1,1) 
mode for gN=731.5 and no = 0.79. Dotted vertical lines 
show frequencies obtained based on the formula 12 U whereas 
dashed ones mark frequencies calculated from the simplified 
expression 1181 . 
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FIG. 7: Dispersion relation for various interaction strengths 
and temperatures. Numerical results are marked by boxes 
(gN = 731.5, n = 0.95), circles (gN = 73.15, n = 0.13), and 
triangles (gN = 73.15, no = 0.93). Solid lines are obtained 
from the expression I21H . Inset shows excellent agreement 
between the numerical data and formula 12 U for all modes in 
the case of high temperature. 



Approximate analytical solutions show that ampli- 
tudes a p and a* oscillate with two frequencies. By tak- 
ing appropriate linear combination of these amplitudes 
one can find normal modes of the system, i.e., modes os- 
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FIG. 8: Quasiparticle energies as a function of inverse of the 
population for the system of N = 235000 (gN = 731.5) atoms 
for various temperatures. The upper set of data (circles) was 
obtained on the grid 36x36x36 and corresponds to the tem- 
perature of 9.37 nK (no = 0.07), the middle one (triangles) is 
the case of 4.27 nK (no = 0.79) and 24x24x24 grid whereas 
the lower set (marked by boxes) results from the 16x16x16 
grid and gives the temperature of 1.80 nK (no = 0.95). 



cillating with only one frequency. This transformation is 
just the transformation to Bogoliubov quasiparticle am- 
plitudes. They can be regarded as an elementary exci- 
tations of the system. However, numerical results show 
that instead of two frequencies one obtains two distinct 
groups of sharp frequencies. For each positive frequency 
component /i + e (where e corresponds to the frequency 
within the peak) there exists companion at the frequency 
(jl— e. Moreover, the phase difference of these components 
is equal to the corresponding one for mode —p. The cen- 
ter of the higher frequencies group defines the energy of 
quasiparticle while the width of this group is related to its 
lifetime. By integrating over the quasiparticle spectrum 
one can get the occupation of the quasiparticle mode. 
In the phonon range the occupation of the quasiparti- 
cle must be obtained with the help of the Bogoliubov 
transformation. This task is simplified because of the 
phase relation mentioned above. Having energies of ele- 
mentary excitations and their populations we can assign 
the temperature to the system 0. Since the modes are 
macroscopically occupied we expect that the equiparti- 
tion relation is fulfilled, i.e., Nn p (e p — fj,) = ksT for each 
mode. As it is the case shows Fig. [H] where we plot- 
ted the energy of quasiparticles as a function of inverted 
population. Due to the equipartition relation this depen- 
dence should be linear and the slope of the line is just 
the temperature of the system. However, since the range 
of the 'y' axis in Fig. [SJis rather large some details for 
small momenta might be hidden. Indeed, the presenta- 
tion of the same results in a different form (see Fig. |5J) 
shows departure of low energy modes (phonons) from the 
equipartition law. Finally, in Fig. HUI we plot the conden- 
sate occupation for various temperatures for the system 



built of N = 235000 atoms. The dashed line in Fig. 
ITTjl shows the population of an ideal condensate of the 
same density. Although we can not precisely determine 
the value of critical temperature of finite size interacting 
gas, Fig. ^| clearly shows that the shift of the critical 
temperature is positive. 



12000 




20 40 60 80 100 120 
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FIG. 9: Verifying the equipartition relation we plot the ex- 
pression Nn p (e p — fi) (in units of h 2 / (mL 2 )) for all modes. 
Note the imperfect equipartition for small momenta (phonon 
modes) . Parameters are the same as in the case of middle set 
of data in Fig. [H] 




T (nK) 

FIG. 10: Condensate occupation as a function of tempera- 
ture for a uniform weakly interacting Bose gas; the scattering 
length equals a s = 25 nm (gN = 731.5, N = 235000). The 
dashed line shows the fraction of atoms in the condensate for 
the ideal gas of the same density (the critical temperature is 
about T = 7.1 nK. The solid line is the fit: 1 - (T/9.85) 1 ' 80 . 

Finally, we would like to discuss the damping of the 
quasiparticle modes. There exists analytical result for 
the decay rates of elementary excitations in a uniform 
Bose gas. At high temperatures (ksT ^> gN/V) Landau 
processes dominate and the decay rate for the phonon- 
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approaches agree quite well except of two first points. 
These points correspond to very low temperatures (ap- 
proximately 2000 and 4000 in units of h 2 /{mL 2 )) and 
certainly the condition of the validity of analytical for- 
mula H'22f> is not fulfilled in this case since gN/V is as 
high as 731.5. However, in our method we can go beyond 
the phonon-like excitations and find the decay rates for 
particle-like part of the spectrum. In Fig. ^]we plot the 
damping rates for all momenta in the case of higher tem- 
perature. The solid line follows the analytical expression 
(|22H . The numerical results and formula (12211 agree only 
for small momenta as it should be expected (see Fig. 1110 . 
The damping rates of particle-like excitations show the 
nonlinear dependence on the momentum. 



FIG. 11: Comparison of the ratio y P /p calculated numer- 
ically (circles) with that coming from the analytical for- 
mula 1221 (and given by boxes) for various temperatures and 
gN = 731.5. 




MOMENTUM (ft/L) 

FIG. 12: The decay rates as a function of momentum for 
gN — 731.5 and higher temperature (no = 0.48). Both 
phonon-like and particle-like parts of spectrum are consid- 
ered. Solid line follows the analytical formula 12211 . 

like part of the spectrum follows the formula 14]: 

,p 2 gN/V \ V V h V ' 

Since phonon energies depend linearly on the momen- 
tum, so do the decay rates. The ratio j p /p is a function 
of temperature given by ^^/n^T, where tiq is the con- 
densate fraction. In Fig. (|11J) we compare numerical 
results with the analytical expression (|22|) . In order to 
determine the damping rates we average the spectrum 
of phonon-like modes over the angles in the momentum 
space and fit to the Lorentzian curve. The FWHM of the 
Lorentzian fit is just the damping rate. We see that both 



V. CONCLUSIONS 

We have applied the classical fields approach to the 
weakly interacting Bose gas in a box with periodic bound- 
ary conditions. By Fourier transform we analyze the fre- 
quency dependence of amplitudes of eigenmodes of the 
system. We observe that the zero-momentum compo- 
nent of the field oscillates with single frequency which 
is identified as a chemical potential. We have obtained 
nonperturbative results for the chemical potential. For 
the physically relevant parameters it differs strongly from 
the simple two-gas formula. All other modes have much 
more complicated time evolution. Their spectra consist 
of two groups of frequencies displaced symmetrically with 
respect to the chemical potential. The center of gravity 
of each groups of frequencies agrees perfectly with a sim- 
ple formula known as the Popov formula 11]. The low 
momenta excitations have a phonon-like dispersion while 
for high momenta components the dispersion is of particle 
kind. Knowing the excitation spectrum of quasiparticles 
we are able to find also their populations and we can test 
the equipartition of the energy between the modes and 
then determine the temperature. The width of the spec- 
trum is identified as the damping rate of quasiparticles. 
We observe that our nonperturbative results agree rea- 
sonably well with the analytic estimates for the phonon 
damping rates for high temperatures. 
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